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The structure of a nonuniform Lennard- Jones (LJ) liquid near a hard wall is approximated by that of 
a reference fluid with repulsive intermolecular forces in a self-consistently determined external mean 
field incorporating the effects of attractive forces. We calculate the reference fluid structure by first 
determining the response to the slowly varying part of the field alone, followed by the response to the 
harshly repulsive part. Both steps can be carried out very accurately, as confirmed by MC simulations, 
and good agreement with the structure of the full LJ fiuid is found. 
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Dense liquids have highly nontrivial density correlations 
arising because the harshly repulsive molecular cores can- 
not overlap 1|-|3). Because of the constantly changing 
molecular arrangements, such correlations play a much 
more fundamental role in liquids than they do in other 
condensed phases such as glasses and solids, which sample 
only a few basic configurations. Indeed, a model with only 
repulsive intermolecular forces Q can give a surprisingly 
accurate description of the full density correlations seen in 
a uniform dense simple liquid like Ar because the vector 
sum of the longer ranged attractive forces on a given par- 
ticle essentially cancels in most typical configurations. 

Nonuniform liquids present a greater and qualitatively 
different challenge, since even the averaged effects of at- 
tractive forces clearly do not cancel We discuss here 
an example where both attractive and repulsive forces can 
greatly influence the liquid's structure: a Lennard-Jones 
(LJ) fluid next to a hard wall. We obtain accurate numer- 
ical results using a physically suggestive generalized mean 
field description of the attractive forces |^ . We consider 
first the effects of these slowly varying forces on the liq- 
uid's structure before taking account of the response to 
the rapidly varying (hard core like) part of the external 
field. This treatment of attractive interactions is quite 
different from that used in conventional integral equation 
and density functional methods |^ , and we believe it offers 
important conceptual and computational advantages. 

Fluid particles interact with a known external (wall) 
field <j){r) and through the LJ pair potential w{rij) = 
uo{rij) + ui{rij), divided as usual ||| so that all the re- 
pulsive intermolecular forces arise from uq and all the 
attractive forces from ui. We assume that the external 
field (^(r) = (j)o{r) + 0i(r) can be separated in a similar 
way, where the subscript denotes in all that follows a 
harshly repulsive interaction and the subscript 1 a much 
more slowly varying interaction usually associated with 
attractive forces. Here we take (/)(r) to be a hard wall 
potential, setting (f>i{r) = and (/>o(r) = (f'Hwiz), where 
(^HW — oo for z < 1 (in reduced units) and otherwise. 



and we let be the bulk density far from the wall. 

We relate the structure of the nonuniform LJ system 
to that of a simpler nonuniform reference fluid [§,0, 
with only repulsive intermolecular pair interactions uo('"y ) 
(equal to the LJ repulsions) in a different effective refer- 
ence field (ERF) (j)R{r). While the replacement of attrac- 
tive pair interactions by an approximate "molecular field" 
is an essential step in mean field theory, we can think of 
other more general choices. Here we determine <^R(r) for- 
mally by the requirement that it has a functional form 
such that the local (singlet) density at every point r in the 
reference fluid equals that of the full LJ fluid B: 



pair; [0r]) = p(r; 



(1) 



The subscript reminds us that the reference system pair 
interactions arise only from uq and the notation indi- 
cates that all distribution functions are functionals of the 
appropriate external field. 

To find (/)rt explicitly, we subtract the first equations of 
the YBG hierarchy for the full and reference systems 
with (j>„ chosen so that Eq. (|l|) is satisfied |^,|^ . The result 
can be written exactly as 



- Vi[0H(ri)- 



(ri)] = - y"rfr2po(r2|ri; [(?!'R])ViUi(ri2) 



dr2{p(r2|ri; [0]) - po(r2|ri; [(f>„])}Viw{ri2) ■ (2) 



Here po(i"2|ri; [0r]) = Po^^(ri,r2; [(/)r])/po(i"i; [^r]) is the 
conditional singlet density, i.e., the density at r2 given 
that a particle is fixed at ri . 

If we assume that Eq. (|l]) produces similar local envi- 
ronments for the (identical) repulsive cores in the two flu- 
ids, which mainly determine density correlations through 
excluded volume effects, then the conditional singlet den- 
sities in the two fluids should also be very similar. This 
key structural assumption introduces a generalized mean 
field theory in which the reference fluid still has nontrivial 
pair and higher order correlations induced by the repul- 
sive forces. This suggests that the last term on the R.H.S. 
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in Eq. (g) is often very small. If we ignore it completely 
j|] we obtain the approximate equation for the field (p^ 
suggested by Weeks, Selinger, and Broughton H: 
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Vi[0R(ri) - 0(ri)] = y dr2po(r2|ri; [0R])ViUi(ri2) 



(3) 



Eq. incorporates mean field ideas, but it appropri- 
ately focuses on forces . The relation to ordinary mean 
field theory becomes clearer 0] if we replace po(r2|ri; [(l>n]) 
by Pa {t2 ; [^n] ) in Eq. (^) . This approximation is much bet- 
ter than one might at first suppose, since the main differ- 
ence in these two functions occurs when V2 is close to ri, 
but then for small ri2 the multiplicative factor — ViUi(ri2) 
(the attractive part of the LJ force) vanishes identically. 
The gradient Vi can then be taken outside the integral 
and Eq. (^) can be integrated. Choosing the constant of 
integration so that the density far from the wall equals ps, 
we obtain the simplified mean field equation |Q : 

0R(ri) - (j){Vl) EE 

(j)s{ri)= I dY2[po{r2\[(t)Ti]) ~ PB\ui{ri2) . (4) 



Because of the integration over the slowly varying at- 
tractive potential "weighting function" Mi(ri2), <^s(r) in 
Eq. (Q) extends smoothly into the repulsive core region 
where po(r; [0r]) vanishes. Outside the wall it is smooth 
and relatively slowly varying even when po{v. [0r]) itself 
has pronounced oscillations. Physically 0s (r) mimics the 
effects of the unbalanced attractive forces in the LJ sys- 
tem, giving a soft repulsive interaction Q that tends to 
push the reference particles away from the wall. 

In order to solve equations like ^ or (^) to obtain the 
self-consistent ERF 4>t^{v), we must determine the required 
reference fluid distribution functions arising from a given 
external field. In previous work computer simula- 

tions were used for this purpose. We now introduce a sim- 
ple and accurate numerical method for calculating these 
distribution functions and illustrate it here by solving 
for the case of the LJ fluid near the hard wall. 

We note that the ERF (/)R(r) = (/)Ro(r) -I- (t>ni{r) in Eq. 
(IJ) (and other related equations) can be naturally sepa- 
rated into the sum of a harshly repulsive part, 0Ro(r), and 
a much more slowly varying "smooth" part 4>^i{r), aris- 
ing physically mainly from the attractive interactions in 
the original system. Eq. (||) suggests the identification 
(/)Ro(r) = 00 (r) and 0Ri(r) = ^^(r) -f 0i(r). More gener- 
ally, we can define 4>no{v) = 4>o{v) ~ 0os(i") and 0Ri(r) = 
(/)s(r) -I- 00s (r) + 01 (r), where 0os(r) is an essentially arbi- 
trary smooth function that is nonzero only in the repulsive 
core region but with 0os (r) < < 0o (r) , so that 0ro remains 
a harshly repulsive interaction. In the present case it is 
sufficient to take the separation suggested by Eq. (^), with 
0ro = <i>HW and 0ri = 0^. 

Our task is now to determine the local density 
Po(r; [0r]) = po,R(r) produced by a given ERF 0r. We 
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FIG. 1. Self-consistent potential 0ri = (ps (dashed line) for 
PB = 0.785 (a) and ps = 0.45 (b), T = 1.35, and bare wall 
potential 0ro = <f>Hw (solid line). The ERF 0r = 0ro -|- 4>ri- 

provide a new way to solve this basic problem, quite in- 
dependent of its origins in the mean field equation (|^). 
Initially we treat the LJ repulsive potential uq as a hard 
core interaction, but then use standard methods ||] to 
correct for its finite softness in our final numerical results. 
We expect that there will be very different responses of 
the reference fluid density to the rapidly and slowly vary- 
ing parts of the ERF 0r = 0ro -I- 0ri and anticipate that 
any large oscillations arise mainly from the harshly repul- 
sive part 0ro. These oscillations cause problems in density 
functional methods, which use a variety of somewhat arbi- 
trary "weighting functions" to arrive at some underlying 
"smooth density" for use in a free energy functional ||^. 

Instead, we first determine the response to the slowly 
varying part of the ERF alone, followed by the response to 
the harshly repulsive part, using different methods in each 
step appropriate for the different density responses. In the 
first step, we determine the associated "smooth interface" 
po(r; [0Ri]) = /Oo,Ri(r) that arises naturally from the slowly 
varying part 0ri — (ps of the ERF alone. Physically, this 
takes account of the effects of the attractive interactions 
modeled by 0ri. We start from the basic equation relating 
small changes in the potential and density 

-/3(S0Ri(ri)= / dr2Xo ^(ri,r2; [po,Ri])'5/3o,Ri(r2) , (5) 



through the generalized linear response function 
Xo ^(ri,r2; [po.Ri]) = <5(ri-r2)/po,Ri(ri)-co(ri, r2; [pcri])- 
Here co(ri,r2; [po,Ri]) is the direct correlation function of 
the reference fluid with density po,Ri and [3 — l/ksT. 
Specializing to the case when the change is a small dis- 
placement of the field, we find the exact equation |10|: 



Vi/9o,Ri(i"i) / /Oo,Ri(i"i) = -/3Vi0Ri(ri) 

-I- y dr2Co(ri,r2; [po,Ri])V2PQ,Ri(r2) . (6) 

If Po,Ri is relatively slowly varying, we can accurately 
approximate co(ri, r2; [po.Ri]) under the integral in Eq. 
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(|^) by the uniform fluid function co(ri2;pi2), where pi2 
is some intermediate density associated with the two 
points A natural choice that gives very good re- 

sults when po.Ri is reasonably smooth is pi2 = [po,Ri(ri) + 
Po,Ri(r2)]/2. Starting with a given (/)ri, we can then solve 
Eq. for the associated po.Ri by iteration, making use 
of the analytic and accurate Percus-Yevick (PY) expres- 
sions for the direct correlation function of the uniform hard 
sphere fluid If necessary, we can choose 0ri in- 

side the repulsive core region to help ensure that po.Ri is 
smooth enough for the expansion method to be accurate; 
this procedure is important in some other applications ||] . 

A special case where this step can be carried out analyti- 
cally arises when po.Ri varies so slowly that it is accurate to 
keep only the first term in the expansion of V2/9o,Ri(r2) in 
Eq. (H) about ri. After integrating, we arrive at the sim- 
ple local hydrostatic relation between po.Ri and (pui- 

Mo(po,Ri(r)) + (/!)Ri(r) = ^qb , (7) 

where fJ.o{p) is the chemical potential of the uniform (hard 
sphere) reference fluid at density p and pqb = Po{pb)- 

The smooth profile po.Ri is analogous to one that could 
be calculated using a single occupancy lattice gas (Ising) 
model, where correlations arise only from attractive inter- 
actions [||. A realistic fluid has additional short wave- 
length correlations due to the repulsive intermolecular 
forces. These show up primarily in the second step of our 
method, where we take account of the response to (/)ro, the 
remaining harshly repulsive part of the ERF. 

Consider flrst a small perturbing potential (5(/iro that is 
nonzero only inside the wall region with z < 1. Evaluating 
Eq. (^) for zi > 1 gives an exact relation between the 
small induced density changes inside and outside the wall 
region. However, it has been shown that even large density 
fluctuations in a hard sphere fluid are accurately described 
by gaussian fluctuation theory This suggests that if 
we could somehow impose the proper values on the wall 
density fleld for z < 1 arising from the full ^rOj we could 
then still use the linear response relation to determine the 
large density change Apo^^{r) = po,R(r) — po,Ri(r) induced 
for z > 1. Imposing accurate density values in general is 
very difficult |l^, but for the hard wall potential (/)ro = 
4>HW we have the exact result po,R(r) = for all z < 1. 
Thus replacing 5/5o,ri by Apo.R in (||) and setting po.R = 
for all z < 1, we find for zi > 1: 

Apo,R(ri)/po,Ri(i"i) = Jdr2 co(ri, r2; [po,Ri])Apo,R(r2). 

(8) 

Eq. (|8|) is a linear equation for Apo,R(ri), which we can 
directly solve by iteration or other means, approximating 
Co(ri, r2; [po,Ri]) by that of an appropriate uniform sys- 
tem, just as we did before. When po,Ri(r) = pn, Eq. 
(H) is equivalent to the usual hard-wall, hard-particle PY 
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FIG. 2. Density profiles po,Ri (dashed line) and po,R (solid 
line) compared to MC simulations of the reference fluid in po- 
tential (j!)r (circles) for the same states as in Fig.l. 

equation, which has an analytic solution [Q. Eq. (||) is 
quite adequate for our purposes here, though small errors 
can be seen at the highest densities. If still more accuracy 
is required, we could use modified GMSA type equations 
related to the PY equation It may also be possible 
to use new and very accurate density functional methods 
for hard core fluids in this step of our method ]l5| ] . 

The net result of this two step process is the desired 
po,R arising from a given 0r. This can be substituted into 
Eq. (^ , which can then be iterated to determine the final 
self-consistent In Fig. (1) we give the self-consistent 
potentials 0ri — (jj^ that satisfy Eqs. (^), (^), and (^) 
for two different states along the near critical isotherm 
T = 1.35. We see that (/)ri is indeed a slowly varying 
repulsive interaction in both cases. In Fig. (2) we give 
the associated smooth density profiles po.Ri from (^) for 
each state, as well as the full profiles po,R determined from 
Eq. (^). These are compared to Monte Carlo simulations 
we carried out Q of the reference system in the ERFs of 
Fig. (1). This directly tests the accuracy of our two step 
procedure for calculating the effects of (f>n on the reference 
system. The agreement is excellent. 

In Fig. (3) we test the simplified mean field treatment 
of the attractive interactions in Eq. (^ by comparing the 
reference profiles po,R to those of the full LJ fluid in the 
presence of the hard wall, as determined by MC calcula- 
tions. There is good agreement, though small quantitative 
differences can be seen. Thus even the simplest mean field 
treatment of attractive interactions is capable of captur- 
ing the major changes in the density profile as the density 
is decreased, and at lower temperature at coexistence we 
find complete drying states where a stable vapor-liquid 
interface can exist arbitrarily far from the wall j^]. 

Our emphasis thus far has been on quantitative numer- 
ical calculations. However, the qualitative features of our 
method are equally important. A long-standing problem 
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FIG. 3. Density profiles po,R compared to MC simulation of 
the full LJ fluid (circles) for the same states as in Fig.l. 

of liquid state theory, well illustrated by the nonuniform 
fluid example studied here, is how to treat consistently the 
oscillating molecular scale "excluded volume" correlations 
arising from repulsive intermolecular forces and the more 
slowly varying and longer ranged correlations arising from 
attractive forces and often associated with the formation 
of interfaces. In principle these issues are addressed by 
modern density functional and integral equation methods, 
but in practice a number of uncontrolled and often math- 
ematically motivated approximation are made. It is often 
difficult to assess their physical implications, and to deter- 
mine where the major sources of error lie. Here, we have 
divided this problem into several distinct parts, whose ac- 
curacy can be examined separately, and where the physi- 
cal content and limitations of the methods used are more 
clear. We used here the simplest mean field equation to 
determine the ERF (pj^ but more accurate (though more 
complicated) equations derived from (||) are available ||] . 

To determine the structure of the reference fluid in the 
presence of a given ERF, we first calculated the response 
to the slowly varying part of the ERF, by expanding 
about a uniform system. Next we used a gaussian field 
model [Q (equivalent to a modified wall-particle PY equa- 
tion) to calculate the response to the remaining harshly 
repulsive part (pj^o of the ERF. More accurate methods 
could be used in both steps if necessary, and for qualita- 
tive purposes both steps can be simplified considerably. 
For example, Lum, Chandler, and Weeks have devel- 
oped very simple approximations for use with continuum 
Landau-Ginsburgh type equations that give good qualita- 
tive results in a number of different cases, including hy- 
drophobic interactions in water. Application of these ideas 
to a variety of different problems is underway ||^,[l6|| . 
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